Numerical studies of domains and bubbles of Langmuir monolayers 
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A numerical algorithm based on the finite element methods has been developed to accurately 
determine the shape of the boundary of a domain containing "boojum" textures. Within the context 
of the simple model we adopt, the effects of both bulk elastic anisotropy and line-tension anisotropy 
on the domain boundary can be examined. It is found that line-tension anisotropy must be present 
in order to account for domains with protruding features. Both elastic anisotropy and anisotropic 
line-tension can result in domains with indentations. The numerical algorithm has been extended 
to investigate the problem of a bubble in extended region ordered phase. 
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I. INTRODUCTION 



A Langmuir monolayer is a single molecular layer of 
insoluble surfactant molecules spread on the air/water 
interface. The surfactants are typically amphiphilic 
molecules with a hydrophilic headgroup and a hydropho- 
bic tail. Each of the individual molecules has internal 
degrees of freedom, namely the tilt and the tilt azimuth. 
Such a system exhibits complex phase structure . The 
"tilted" phases have uniform tilt and possess mesoscopic 
ordering in the tilt azimuth. The structure of the tilt 
azimuth is typically observed as a variation in the light 
intensity under a polarized light microscope. The tilt az- 
imuth organization is referred to as the texture. Various 
classes of the texture have been observed, such as the 
stripes in the bulk ||], the star configuration |^,^ and 
"boojums" Q in the domains of the tilted phase, when it 
coexists with an isotropic phase. The term Boojum refers 
to a class of textures that has a tilt azimuth distribution 
which resembles the structure of the orbital angular mo- 
mentum in a superfluid ^He droplet [^. The "boojum" 
texture, in which the tilt azimuth is distributed continu- 
ously without singularity, will be the subject of this re- 
port. Domains observed to contain a boojum texture are 
not circular in shape In addition, micron-size bub- 

bles, which are regions of isotropic phase surrounded by 
a "tilted" phase, have been found to have non-circular 
shapes ||^. The local tilt azimuth in the "tilted" phase 
around the bubble exhibits a non-trivial structure, which 
has been termed an "inverse boojum" . The relationship 
between experimentally observed textures and the under- 
lying structure of the ordered phase has attracted atten- 
tion in the literature recently. In particular, the "boo- 
jum" texture was first discussed by Mermin in the con- 
text of orbital angular momental distribution in super- 
fluid ■^He droplet Similar texture of has been found 
and discussed in the liquid crystal films l^],^. An exten- 
sive discussion on the various classes of textures in the 
Langmuir monolayers can be found in Ref. 

The problem of the equilibrium shape of, and the tex- 



ture contained in domains in a Langmuir monolayer has 
been investigated by Rudnick and Bruinsma who 
varied both the texture and the boundary analytically 
in a perturbative manner. It was discovered that a non- 
circular boundary represents the equilibrium shape of a 
domain only when there is bulk or line-tension anisotropy 
The equilibrium domain boundary was derived 
as a function of line-tension anisotropy. Galatola and 
Fournier obtained numerically, in a fixed background 
texture, the equilibrium boundary when both elastic and 
line-tension anisotropics are present. Riviere and Me- 
unier have attempted to explain the observed non- 
circular domains in terms of elastic anisotropy. Evidence 
of bubbles with a non-circular boundary and an "inverse 
boojum" has been reported, and a qualitative theoretical 
discussion of the equilibrium shape and texture associ- 
ated with the bubbles can be found in Ref. 0. In the 
spirit of Ref. jl^ , the authors have analyzed in Ref. |T^ ] 
the equilibrium texture and boundary shape combina- 
tions perturbatively to the first order in both the bulk 
elastic and line-tension anisotropics. The approach de- 
scribes the infinitesimal response of the texture and the 
boundary to anisotropic parameters. However, when the 
correction is large enough to be observed, the validity of 
first order perturbative calculations becomes question- 
able. The extension of the perturbative approach to in- 
clude higher order corrections is algebraically formidable. 
If one is to go beyond first order effects, the use of nu- 
merical techniques in this problem is inevitable. 

The major challenges in this problem are, first, the 
evaluation of a 2D texture with a boundary condition 
on the boundary, which is, itself, variable. Secondly, not 
only must the texture be evaluated with high accuracy, 
but a precise determinations of the derivatives of the tex- 
ture on the boundary is also crucial to the computation 
of the boundary shape. The authors have developed a 
numerical algorithm based on the finite element method 
(FEM) with adaptive mesh refinement for the evalu- 
ation of a 2D texture and its derivatives. The boundary 
corrections can then be computed using the Runge-Kutta 
methods llj] . Implementation of the numerical method 
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reveals various classes of domain shapes ranging from 
those with indentations to those with protruding features 
and, additionally, of cigar-shaped domains. The effects of 
bulk elastic anisotropy have also been examined. These 
studies lead us to the conclusion that, as least for those 
domain shapes observed to date, it is more likely that 
the line-tension anisotropy is responsible for non-circular 
domains. A brief account of the study described above 
has appeared in an earlier publication . The numeri- 
cal results also confirm that the qualitative conclusion to 
be drawn from the perturbative treatment are preserved 
up to large anisotropic parameter. 

In this report, we describe in detail the implementa- 
tion of the numerical methods that lead us to the results 
reported in Ref The extension of the algorithm to 
allow for computation in the case of bubble has also been 
examined. It is verified numerically that bubbles acquire 
a non-trivial boundary shape when only the first term in 
the Fourier expansion of the line tension is present. This 
result contrasts with what is known to be true in the 
case of domain, which remains circular in the presence 
of this low-order line-tension anisotropy With the 

use of our numerical algorithm, we are able to examine 
the effects of the bulk elastic anisotropy on the shape 
of the bubble and on the texture that surrounds it. We 
find that bulk elastic anisotropy significantly affects the 
texture in the condensed phase around the bubble while 
leaving the boundary nearly unmodified. 

The organization of this paper is as follows. Section || 
contains the details of the computational scheme for the 
evaluation of the equilibrium textural and boundary con- 
figuration for domains. The discussion covers the deriva- 
tion of the simplest variational formulation of the finite 
element method in our specific application, the Runge- 
Kutta method and the combined algorithm. In Sec. [11, 
results for the domain are examined. Section [V describes 



the extension of the numerical algorithm to the problem 
of bubbles. An examination of the results of the pertur- 
bative treatment follows. New results on the effect of the 
bulk elastic anisotropy on the textures around the bub- 
bles are discussed. Finally, Sec. ^ contains concluding 
remarks and discusses possible future extensions of the 
numerical methods discussed in this report. 



II. THE NUMERICAL ALGORITHM 



The model that we adopt for the Langmuir mono- 
layer is a simple elastic model of an ordered media 
associated with XF-like order parameter — a 2 dimen- 
sional unit vector c{x, y), which can be parameterized as 
xcosQ{x,y) + ysiiiQ{x,y) The quantities x and 

y are unit vectors in a Cartesian coordinate system, and 
0(a;, y) is the angle between c and the x-axis. The energy 
of the system contains contributions from the boundary, 
r, in addition to the bulk, il. The most general form of 
the elastic energy for such a system with in-plane 



reflection symmetry (an achiral system) can be written 
as 

H[e]= [ nbdA+ (f cr{d-&)ds, (2.1) 
Jn Jr 



where 



2 ' ' 2 
(7((/>) — ao + a„ cos ncj) 



(2.2) 
(2.3) 



Ks, Kb are respectively the splay and bend elastic mod- 
uli, and t? is the angle between the outward normal n of F 
and the x-axis. The setup of the computations is shown 
in Fig. 1^. In terms of the average Frank modulus k and 
the coefficient of elastic anisotropy, b where 2k = K^ + Ks 
and 2k6 = Ks — A'b, the extrema of the elastic energy 





FIG. 1. The geometry of the calculations for (a) domains 
and (b) bubbles in plane-polar coordinates where the bound- 
ary r is parameterized by pif). The gray area is the bulk 
designated by Q. n and t are the outward normal and the 
tangent, respectively. Q is the angle between the c-director 
and the x-axis and i9 is the angle between the outward normal 
of the boundary and the x-axis. 
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Eq. ( ^.l| ) occurs when Q{x,y) and the bounding curve T 
satisfy their respective equiUbrium conditions. The ex- 
tremum equations for Q{x,y) are 

- V^e + b [(e^^ - Oyy) cos 29 + 26^,^ sin 26 

+ (-e^ + e^) sin 29 + 29^9y cos 29] = 0, (2.4) 

in Q and 



k9„[1-6cos2(?9-9)] + 
KbQt sin 2(i? - 9) - cr'(i9 - 9)= 0, 



(2.5) 



along r, where 9„ = n- V9, 9i = t- V9, t being the tan- 
gential vector. The extremum equation for the bounding 
curve r, in terms of 9„, 9f and (I'd/els, is 



Hb~<j'{^-e)en-a"{'0 

+ [cr(i?-9) + a"(i? 



9)9* 



(2,6) 



where ds is the length element of F traversing in the pos- 
itive direction of 17 and A is a Lagrange multiplier that 
enforces the condition of constant enclosed area. 

The equations for both 9 and F are complex and 
highly non-linear. Closed form analytic solution of the 
extremum equations is almost impossible. Attempts have 
been made to solve the simultaneous equation pertur- 
batively to first order in the elastic and line-tension 
anistropies |12|. When the corrections to the bound- 
aries are large enough to be observable, it is not expected 
that the results are accurate and high order corrections 
have to be taken into account. However, these attempts 
provide us with insight with regard to the infinitesimal 
response of the boundary to the anisotropics under inves- 
tigation. In the work to be described below, we analyze 
the equations numerically in order to further explore the 
implications of the simple model Eq. (2.1) for a larger 
range of the anisotropic parameters. We retain coeffi- 
cients up to 02 in the expansion of the line tension in 
our analysis, i.e. <j{(t)) = (Tq -I- ai(cos -I- 7 cos 20), where 
the quantity 7 = a-ijax is defined for convenience. We 
remark that the analysis will be based on the exact "boo- 
jum" texture with circular domain when 7 = 6 = 0. The 
boundary will be computed in terms of the corrections to 
the circular boundary. The discussions will be restricted 
to those domains with boundaries F for which the dis- 
tance from each points on the curves to the origin e'^^'^^ 
is a single- valued function of the polar angle ip. 

The numerical algorithm consists of two parts: in the 
first part, one evaluates the texture 9 using an assumed 
boundary F, and, in the second part, one computes F 
using a fixed 9. Simultaneous equilibrium conditions 
for F and 9 are achieved when a set of predefined self- 
consistent criteria is met. It is evident from the form 
of Eq. (|]^) that accurate determinations of 9 and its 
derivatives are the key factors in the solution of the prob- 
lem. The requirement that the assumed F in the first part 
of the algorithm be an arbitrary curve rules out finite dif- 
ference methods, and militates in favor of finite element 



methods (FEM). A key feature of the FEM is flexibility 
in the choice the set of points at which the functional val- 
ues are to be evaluated, including those on the boundary 
of the region of interest. This feature is exactly what is 
needed in our problem, because of the non-trivial geome- 
try of the boundary. Cne of the simplest constructions of 
the FEM in 2 dimension is described as follows jl^ . We 
first approximate F by a polygonal curve, then subdivide 
n into a set of non-overlapping triangles. No vertex of 
one triangle lies on the edge of another in the set. The 
edges of the set of triangle forms a mesh that covers fi. 
The process of creating this set of triangles is called mesh 
generation. The resulting set of triangles is referred to 
as the triangulation of 17. Functions are defined by their 
values on the vertices of the triangles in the triangulation. 
The value of a function within a triangle is obtained by 
interpolation using the values on the vertices. Integra- 
tion over r2 is the sum of integrations over the triangles 
which can generally be trivially evaluated. We have now 
projected our problem, originally on an infinite dimen- 
sional space onto a N dimensional space, where N is the 
number of vertices in the triangulation of f2. We may 
write = (9i), %= 1, • • • , N and 



N 



(2.7) 



where ipi{x, y) is a set of basis functions of the N dimen- 
sional space. These ipi 's should not be confused with the 
polar angle which is denoted by the symbol without 
a subscript. The discrete version of the elastic energy 
functional Eq. (2.1) is a function of N variables 9^ and 
it can be rewritten in terms of k and h as 

^(0) = f {|V9p + h [{-Ql + 9^) cos 29- 

2Q^Qy sin 29]} dA+ <j> cr{d - 9) ds, (2.8) 

where Ox = J2 ^iVix, = ^i^iy^ Vix = dipt/dx and 
Lpiy = d^pi/dy. The equilibrium condition becomes 



dH{@) 
99,: 



0, 



1, 



■,N, 



(2.9) 



which is a discretized version of Eqs. (2.4) and (2.5). The 
set of above equations is not linear. However, if we write 
them in the form of A(0) • = b(0) where A(0) is an 
N X N matrix, b(0) and are 1 x N column matrix as 
shown below, 

Aij{@) = K / [pixP'jxi'^ - 6 cos 29) 4- 
ipiyPjy{l + 6 COS 29) 4- 

K'Pix'P'jy + Piy^P]x)s,V[i2Q]dA (2.10) 



h{@) 



Kb I [{-ei 



9^) sin 29 
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2Q^Qy cos2Q](p.,dA 
a\^-e)ip^ds 



(2.11) 



we are able to solve for iteratively using a standard 
numerical algorithm for the solution of systems of linear 
equations. We have adopted the method of LU decom- 
position for solving 0. 

The mesh generation algorithm plays an important role 
in the efficiency of the FEM. An adaptive mesh genera- 
tion algorithm is used in our program to determine 0. 
We start with a mesh that is nearly regular throughout 
n with a predefined grid size. After obtaining a first esti- 
mate of 0, a refined mesh is generated. The refined mesh 
has variable grid sizes over f2 depending on the variation 
of 0. Figure || depicts the process of mesh generation 
with adaptive refinement. We are able to determine not 
only 0, but also the derivatives Qt and 8„, which are 
necessary for the evaluation the bounding curve, accu- 
rately and efficiently with the adaptive mesh generation 
algorithm. 

The next part of the algorithm is the determination of 
the bounding curve F. We assume the order parameter 
field 8 is fixed in Eq. ( |2.6| ) so as to simplify the prob- 
lem. We then pick an origin in and parameterize the 
bounding curve F as k{ip), where e''^'^'' = |r| is the dis- 
tance between the origin and {x,y) on F, and ip is the 
polar angle. In this parameterization, Eq. ( ^.6| ) is a sec- 
ond order non-linear differential equation in k{ip). If we 
rewrite Eq. ( p.6D as 

(2.12) 



k" + q{tp;k,k')k' ^r{tp;k,k') 



where k' = dk/dip and 



q{ip;k,k) = ^ — [1 + k ), 



r{ip;k,k') 



a + (7 
a'Qk + fT"8 



(2.13) 



a + a" 



(l + fc'2) + 



e'^Hb + A 



a + cr" 



(l + fc'2)3/2. 



(2.14) 



Again, it is possible to integrate the equation for k(ip) 
iteratively using standard method for the solution of or- 
dinary differential equation. The Runge-Kutta method 
Q is chosen for our application. 

The problem of solving Eqs. (|2^, (O) and (|2^) 



for Q[x,y) and F is refor mulat ed in terms of the so- 
lution of Eqs. ( ^.9[) and ( 2.12 ) iteratively for and 
F. We begin by assuming an initial boundary F'"' 
and texture @^^\ from which the texture 0'^-'^' can be 
computed using FEM. Then iterated texture 0^^^ = 
0(0) + (0'(i) _ 0(o))/i,^o) ^^^^ ^gg^ evaluate a 

new accepted boundary F^^) = F^o) + (F'(i) - F(*")/4°\ 
where F'^^' is obtained using the Runge- Kut ta ordinary 
differential equation integrator on Eq. ( p.6[ ). The pro- 
cess is repeated until both A0(") = |0("nr0("-i) | and 
AF(") = |f(") -F("~i) I are less then a preset tolerance of 
the order O(10~^). The factors u!^^ and with initial 
magnitude in the order of 0(10) is introduced to avoid 
numerical instability in these iterative processes. These 



factors, v, 



in) 



and are adjusted depending on A0(") 
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FIG. 2. The progressive steps towards mesh generation of the problem for the domain, (a) initial choice of grid points in a 
square lattice, (b) deforming to grid points into the region of interest Q and the bounding curve F while keeping the square 
lattice topology and connectivity, (c) triangulation using the square lattice connectivity and (d) resulting mesh after adaptive 
mesh refinement. 



4 



and Ar*^"), and at the final iterations take on values close 
to unity. 



III. DOMAINS 



It is well established that the boundary is strictly cir- 
cular for a domain with a boojum texture when elastic 
anisotropy and line-tension anisotropy are not present, 
or 7 = 6 = JlO| . This texture-bo und ary combination 
is indeed a local minimum of Eq. {2A_) [Q. Thus, in 
order for there to be non-circular domains, it is neces- 



sary to retain terms in the expansion Eq. (2_^) up to 
at least terms going as 02. Using the numerical algo- 
rithm described above, we have performed systematic 
studies of the domain textures and shapes in terms of 
the elastic anisotropy, the line-tension anisotropy as well 
as the domain size. Before we describe our observation, 
we note that when 7 = 5 = 0, the exact result is given 
by a circular boundary of radius Rq together with a 
"boojum" texture with a +2 defect located a distance 

Rb = Rq (^1 + \/l + Po) /pa from the center of the do- 
main, where po = Roai/n is the normalized domain ra- 
dius | [Tc[ | . An exact equilibrium texture-boundary combi- 
nation is shown in Fig. ^. The simulated image obtained 
using the Brewster angle microscopy (BAM) is also dis- 
played in the background. The signature of a "boojum" 
in a BAM image is a set of straight constant-intensity 
lines emerging from a virtual defect slightly outside the 
domain. The light intensity in a BAM image depends on 
the exact experimental setup and the properties of the 
monolayer In the case of all simulated BAM im- 

ages presented in Fig. ^ and elsewhere in this report, the 
Brewster angle is taken to be that of water — 53.12°, 
the angle of the analyzer a is equal to 90° , the thickness 
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FIG. 3. Circular domain with a "boojum" texture, iio = 5, 
ctq = 4, ai = 1.6, 02 = 0, K = 1 and 6 = 0. Also shown on the 
background is the simulated image that would be obtained by 
Brewster angle microscopy. 



of the monolayer is assumed to be d = 0.3nm, the tilt 
^I^ is 30°, the dielectric constants of the monolayer are 
€\_ = 2.31, e|| = 2.53 and it is assumed that the wave- 
length of the light A = 514nm. 

We first concentrate on the effects of h and keep 7 = 0. 
When 6 < 0, the texture is altered in such a way that 
the virtual defect appears to move closer to the bound- 
ary. This is observed as accelerated convergence of the 
constant-intensity lines to a point on the boundary. On 
the other hand, when & > 0, the texture relaxes as if the 
virtual defect has moved away from the boundary. The 
deviation of the texture from the boojum texture is as 
large as 20% when |&| « 0.8. The textural response is 
qualitatively in accord with that reported in Ref. [ pT| . 
Although there are significant textural corrections due 
to the presence of bulk elastic anisotropy, the resultant 
textures very much resemble a "boojum" as seen in a 
BAM image as shown in Fig. IJ. This means that it is 
difficult to identify elastic anisotropy based on observa- 
tion of the textures in the domains. The response of the 
boundary to elastic anisotropy that we obtain contrasts 
to that reported in Ref. |ll| . The domain acquires an 
indentation when 6 < 0. The indentation remains ob- 
servable for a large range of domain sizes. The boundary 
protrudes slightly for & > 0, as depicted in Fig. ^. The 
protrusion when h = +0.8 is subtle and does not resemble 
the sharp feature observed experimentally. Thus, elastic 
anisotropy alone is not capable of accounting for the 




FIG. 4. The textures and the shapes of domains with 
-Ro = 5, (To = 4, a\ = 1.6 and ai = 0. Their stiffness co- 
efficients are (a)/3 = -0.8, (b)/3 = 0.8. 




FIG. 5. The domain shapes computed for ao = 4, ai = 1.6 
and a2 = 0, and i?o = 0.5, 1, 2, 4, 8. Their stiffness coeffi- 
cients are (a)/3 = —0.8, (b)/9 = 0.8. For ease of observation, 
domains are not shown to scale. 
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FIG. 6. The textures and shapes of domains computed for 
-0.5 < a2 < 0.5, Ro = 5, cro = 4, ai = 1. (a)a2 = -0.5, 
(b)a2 = -0.3, (c) a2 = -0.1, ( d) a2 = 0.1, (e) 02 = 0.3 and 
(f) a2 = 0.5. 

shapes of the domains observed experimentally. Figure |^ 
shows domains of various sizes when (a) b — —0.8 and 
(b) b = +0.8. 

We now proceed to discuss the role of the line-tension 
anisotropy, parameterized by 7, in the textures and the 
boundaries of the domains. Elastic anisotropy will be 
eliminated (5 — 0) for simplicity. We first investigate 
situations when I7I < 1. For very small domains where 
Roai/K <^ 1, the texture is almost constant and the dom- 
inant contribution to the boundary deformation comes 
from the a2 contribution. The domain is elongated at 
both ends along the axis connecting its center and the 
virtual defect when 7 > and is flattened at both ends 
along the same axis when 7 < 0. Domain shapes ex- 
hibit a 2-fold symmetry. When RqQi/k > 1, the texture 
closely resembles the boojum texture and contributes sig- 
nificantly to the boundary distortion through the influ- 
ence of 7 in the line-tension. The domain nolonger dis- 
plays 2-fold symmetry and acquires a protrusion when 
7 > 0, or an indentation when 7 < 0. Figure ^ shows 
domains of with 7 ranging from —0.5 to 0.5. The numer- 
ical algorithm also allows us to examine domain shape 
and texture when ai = and 02 = 1. In this case, the 
domain acquires a "cigar-shape" and the texture is asso- 
ciated with two virtual -fl defects ||l^,|l^. The progres- 
sive changes of the texture and the shape from a system 
with oi = 1 and a2 = to one for which ai = and 
02 = 1 are shown in Fig. ^. When both ai and 02 are 



FIG. 7. The textures and the shapes of domains with p — 1 
and ao = 4. Their coefficients of the anisotropy line tensions 
are (a)ai — 1, a2 = 0, (b)ai — 0.8, 02 ~ 0.2,(c)ai = 0.6, 
a2 = 0.4,(d)ai = 0.4, 02 = 0.6,(e)ai = 0.2, 02 = 0.8 and 
(f)ai = 0, a2 = 1. 

nonzero, the texture can be thought of a superposition of 
pure fli and pure a2 textures. Typically at RqUi/k ^ 1, 
the effect of the set of two -1-1 defects become observable 
when 7 = 1/4. Domains with indentation, protrusions, 
and the "cigar-shaped" domains, have all been observed 
experimentally p^ ]. 

We have already briefly discussed the issue of size de- 
pendence in the previous paragraph. To look into this 
matter in detail, we will examine the particular set of 
data reported in Ref. in which the domains investi- 
gated possess protruding features sharp enough that "ex- 
cluded angles", Vl/o, characterizing the boundaries can be 
identified. The definition of '^q and the experimental 
data are depicted in Fig. ^. The key features of this set 
of experimental data are: (1) goes through a maxi- 
mum as Rq varies; (2) there is an abrupt onset of ^'o in 
the small Rq region; and (3) the intercept at the 4'o-axis 
when the curve is extrapolated implies lim ^ 0. 

-Ro^oo 

Before we make comparisons between theoretical re- 
sults and the experimental data, we comment on the ex- 
traction of V&o from computed domain boundaries. It has 
been shown that within the parameter regime of our dis- 
cussions, the domain boundaries are smooth and contin- 
uous. There is no cusp-like singularity on the boundary. 
This can be seen in the domains of various sizes shown 
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in Fig. ^. Nevertheless, ^'o can be unambiguously mea- 
sured for some of these domains. The values of the pa- 
rameters utilized here are k = 1, S ~ 0.4 and 7 — 0.5. To 
determine for these domains, we adopt a systematic 
scheme that utilizes the function / = /oCxp[— (d^x/dy^)^] 
to capture the most likely ^0 for a given domain bound- 
ing curve devised in Ref. [ p^ , where x{y) parameterizes 
r in Cartesian coordinates system. Density plots of / as 
a function of = — 2tan~^ dx/dy and Rq for numerical 
and the perturbative results are shown, respectively, in 
Figs. ^0|(a) and (c) , the darker regions representing larger 
/, and highlighting the more likely values of ^'q. With 
the use of this method for the determination of \l/o, we 
have obtained reasonable agreement between the pertur- 
bative analysis and the numerical computations in the 
large-i?o regime. We note here that the value at which 
7 is set, 0.5, is too large for perturbative results to be 
dependable. However, the perturbative results resemble 
those obtained numerically in the sense that increases 
as Rq decreses from 00. The abrupt onset of 4'o indicated 
in Figs. |l^(c) and (d) is not present in Figs. |l^(a) and 
(b). It is, however, evident in Figs. ^ that can be 
unambiguously identified for domains with -Ro > 1 while 
the domains become elliptical for which \l/o = when 
i?o < 1- Hence, there is an apparent jump in '^q near 
Rq = 1 and beyond which ^'o becomes non-zero. The 



Rod^m) 

1.3 1 Q . i 




l/Ro(l/|^m) 

FIG. 8. (a) The definition of the excluded angle *o- (b) 
Experimental measurements of the domain-size Ro depen- 
dence of 3'o observed in L2 domains surrounded by LE phase 
taken from Ref. [9]. 




FIG. 9. The shapes of domains of various sizes computed 
for (To = 4, ai — 1.6 and 02 — 0.8. (a)Smaller domains with 
Ro = 0.2, 0.25, 0.33, 0.5, 1 which exhibit 2-foId symmetry. 
(b)Larger domains with Ro — 2, 2.5, 3.3, 5, 10 which have a 
protrusion on one end of the boundary. Each of the sets of 
domains are plotted to scale. 



jump in ^'o predicted in the perturbative analysis [T^ is 
indeed confirmed by the more reliable numerical compu- 
tations reported here. For very small domains {Rq < 1), 
the shapes are predicted to be elliptical by our numeri- 
cal analysis, in contrast to the prediction of nearly cir- 
cular domains that results from the perturbative anal- 
ysis. The magnitude of 7 that results in breakdown of 
the first order perturbative analysis is the key origin of 
the mismatch. In figures 0(b) and (d), the maximum 
Imax of / is shown as the dark line segments and the 
grey bands mark the regions in which / > I^axl'^- They 
depict, respectively, numerical and perturbative results. 
Superimposed are the experimental data which provides 
a reference for the comparisons described above. 

To compare the theoretical results to the experimen- 
tal data, a length scale is required. The length scale is 
set by the assignment n/ai = 4/iTO when the compar- 
isons is made between the perturbative results and the 
experimental data ||l^ . Except for examining the results 
of more reliable computations, there is no attempt to fit 
the experimental data in this report for reasons to be 
discussed below. We first adopt the same set of param- 
eters, with which the perturbative analysis fits the data 
in the large- i?o regime, for the comparison. It is obvious 
from Fig. |l^(b) that even at 7 = 0.5, the theoretical 



(a) (b) Rg (i^m) 
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(c) (d) Rg (j^m) 
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FIG. 10. (a) Density plot of J as a function of ^ and Ro 
of the numerical results for the domains, (b) Plot of Imax 
and the region in which / > Imax/ 2 as a function of ^ and 
-Ro of the numerical results for the domains. Superimposed 
are the experimental data shown in Fig. H(b) with param- 
eters «/ai = 4/xm, 5 = 0.4 and 7 = 0.5. (c) and (d) are 
corresponding plots for the perturbative results. 
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prediction for the maximum of ^'o is much smaher than 
that observed experimentally. This superficial compar- 
ison between the maxima of V^o implies that 7 is very 
much larger in the system investigated. Detailed com- 
parisons do show excellent agreement for domains larger 
than lOfim. Experimentally observed domains with max- 
imum and small circular domains are not reproduced 
numerically. Attempts have been made to investigate the 
combined effect of elastic anisotropy b and 7. However, 
for 7 of such a magnitude, contributions from the b do 
not affect the qualitative behaviors discussed in this con- 
text. It is thus concluded that although the simple elastic 
modelis not capable of fully addressing the issue of the 
domain size dependence of the shapes, it has successfully 
produced the qualitative features in the ^0 versus Rq plot 
and many nontrivial domain shapes observed in various 
experiments jis) . 

In the large- i?o regime {Rq ^ 1), the boundary cor- 
rections are confined in a small portion of the boundary 
and the domains become nearly circular. Because of the 
rapid texture variations in the immediate vicinity of the 
boundary, associated with the approach to the boundary 
of the virtual defect, we are unable to perform depend- 
able numerical investigations of extremely large domains. 
This leaves open the question of the asymptotic behavior 
of 'So is the Rq ^ 00 limit. 

With the numerical scheme for evaluating simultane- 
ousl y and T. We are able to explore the simple model 
) in a much wider range of the parameter space 



.1 



Eq. (|2 
with 

confidence. Not only does the model account for the 
domains with various features observed in experiments, 
it also yields an appropriate domain size dependence of 
the boundary shapes. However, we are unable to perform 
reliable numerical investigations on extremely large do- 
mains. Despite the fact that there is an upper bound to 
the domain size that we are able to compute, we believe, 
on the basis of measurements of the defect positions ||] 
that the largest domains we are able to compute are not 
smaller than those that have been observed experimen- 
tally. The numerical algorithm appears to be capable of 
evaluating domain shapes for arbitrary anisotropic line - 
tensi on, w ith one caveat. A closer look at Eqs. (2.13) 
and (2.14) immediately indicates that this approach is 




not appropriate for situations in which a + a" = a,t 
some points on the boundary. An approach that is ap- 
propriate to this situation is the Wulff constructions pOl . 



IV. BUBBLES 



We now turn to the investigation of bubbles. The first 
task is to numerically evaluate the texture in a region, 17, 
that does not have an external boundary. It is possible to 
implement a straightforward extension of the problem of 
the domain by introducing an artificial external bound- 
ary far away from the inner bounding curve F. One must 
introduce a boundary condition on this added external 
boundary by hand. Figure pi] displays the triangulations 
associated with such implementation of the approach to 
the calculation of the property of a bubble. The method, 
though inefficient, produces results that are consistent 
with those obtained perturbatively |]l2| . 

The problem that involves an infinite fl with internal 
boundary F is referred to as the exterior problem. If one 
does not introduce an artificial external boundary, it is 
necessary to have at hand a complete set of exterior so- 






FIG. 11. Mesh with additional external boundary used in 
the direct extension of the FEM algorithm for the bubbles. 



FIG. 12. Solving the problem of the bubble by transform- 
ing it into an inverse domain problem, (a) shows the region 
f2 of ordered phase and the boundary F for the bubble, (b) 
shows the transformed region fi' and boundary F'. (c) shows 
the mesh that covers fl' . 
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FIG. 13. Boundaries and textures around the bubbles for 
-Ro = 5, ai = 1, K = 1 and fe = 0. 

lutions to construct the boundary condition at F. For our 
particular case, this is possible when 6 = 0, in which case 
the bulk extremum equation reduces to Laplace's equa- 
tion. The examination of the problem with nonzero 6 is a 
major goal of this investigation, and we are not aware of 
the existence of an appropriate set of external solutions 
in this case. Noting that the order parameter tends to 
a fixed value (0 = for our case here) as r — s- oo, it is 
possible to approach the problem of the bubble using a 
different set of polar co-ordinates, i.e. {r' ^Lp) = (l/r^ip), 
that transform the bubble into a domain of area fi' and 
bounding curve F', shown in Fig. ^ with the following 
"elastic energy" , 



n = ^J {el + el + p [{-el + el) cos2(e - 2^) 

(t(i9 - 9) 



-2e^eySin2(e-2^)]} dA 



R- 



ds 



With the problem transformed, the meshing algorithm 
used for the domain can be applied immediately. We are 
then able to proceed with the investigation of the bub- 
bles with the same efficiency and same accuracy as the 
studies of the domains. 

In the numerical studies that we have performed with 
the use of the transformation above, the results for the 
bubbles reported in the perturbative analysis that 
the boundaries are not circular even when 5 = and 
02 = 0, have been confirmed. Figure |l3| shows the texture 
and the boundary of a typical bubble. In the background 
simulated BAM image, one notes the circular constant- 
intensity lines which identify the "inverse boojum" . The 
numerical algorithm further enables us to obtained equi- 
librium bubble boundaries and textures around them 
when elastic anisotropics are present. As can be seen in 
Fig. the elastic anisotropy leaves the boundaries sub- 
stantially unaffected while significantly changing the ap- 
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FIG. 14. Boundaries and textures around the bubbles for 
i?o = 5, ai = 1, K = 1, (a) with 6 = -0.8 and (b) h = 0.8. 




FIG. 15. Boundaries of bubbles computed for k = 0.16, 
CTo = 1, ffli = 0.16 and (a) i?o = 1, 2, 3, 4, 5, 6, and (b) 
iio = 8, 20, 40. 

pearance of the textures around the bubbles. The BAM 
images are also shown in the same figures. In contrast to 
the case in which 6 = 0, the constant-intensity lines be- 
come elongated perpendicular to the axis connecting the 
center of the bubble and the position of the virtual defect 
when 6 < 0. These lines are elongated in the direction of 
the axis when 6 > 0, as shown. This allows for the deter- 
mination of the sign of h in the Langmuir monolayer by 
examining the BAM images of the bubbles. 

In Figs. |l^, we display the size dependence of the bub- 
ble boundaries. Bubbles appear to be circular when they 
are small (i?o < !)• For large enough bubbles (i?o > 1), 
an "excluded angle" defined in Fig. ||(a) can be iden- 
tified. An approach similar to the analysis of the size 
dependence of the boundary in the case of domains can 
be applied. Figures |l|(a) and (c) compare the density 
plots of / as a function of ^ versus i?o for the numer- 
ical and the perturbative results. The two plots are in 
excellent agreement. When 7 = 0, the texture docs not 
differ significantly from that of the inverse boojum even 
though the bubble is not exactly circular. This contrasts 
to what is seen in domains when 7 ^ 0, in which case the 
textures can deviate significantly from the boojum tex- 
ture. The i?o dependence of ^'o features that are quali- 
tatively similar to those seen in the case of the domains, 
i.e. a maximum and an onset. These features have also 
been observed experimentally |^]. Experimental data is 
shown together with the numerical and perturbative re- 
sults in Figs. |l^(b) and (d), respectively, and all the 
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FIG. 16. (a) Plot of / as a function of 'if and Ro of the 
numerical results for the bubbles, (b) Plot of Imax and the 
region in which / > Imax/ 2 as a function of 4' and Ro of the 
numerical results for the bubbles. Superimposed are the ex- 
perimental observations of gaseous bubbles in L2 phase. The 
experimental data have appeared in Ref. [9]. The parameters 
for the by-eye fit are «/ai = OAfim and 5 = 0.16. (c) and (d) 
are the corresponding plots for the perturbative results. 

results match reasonably well. The by-eye fit has been 
obtained in Ref. |^ and no further adjustment of the pa- 
rameters is made in this investigation. 

We have thus devised a numerical method to approach 
the problem of the bubble that can be implemented with 
the same efficiency as in the problem of the domain. 
There is good agreement between perturbative and nu- 
merical results. We are able to investigate the effect of 
the elastic anisotropy, and our results point to a possible 
means for the determination of the relative strength of 
Kg and K^. 



V. CONCLUSIONS 

We have discussed in this report the implementation 
of a numerical method that leads us to the solution of 
simultaneous equilibrium conditions for the textures and 
the bounding curve of the domains. Using this numeri- 
cal algorithm, we have investigated the influences on the 
textures and the domain shapes of the line-tension and 
elastic anisotropics. Our analysis of this simple model 



reveals that elastic anisotropy does, indeed, result in in- 
teresting domain boundaries with protrusions and inden- 
tations. The domains with indentations resemble those 
observed experimentally. However, the domains that we 
generate with protrusions are very different from those in 
observed in BAM images [|l8|. Hence, elastic anisotropy 
cannot qualitatively account for all experimental obser- 
vations. Furthermore, our numerical results are in con- 
trast to the claims in Ref. |jll|]. Dents in boundaries are 
due to a bend modulus that exceeds the splay modulus, 
i.e. b = [Kg — Kb) /{Kg + Kb) < 0, while protrusions 
are present when 6 > 0. On the other hand, the second 
harmonic contribution to the line-tension, parameterized 
by 7 = a2/ai is capable of producing nontrivial domain 
shapes that resemble the shapes observed experimentally. 
For the influence of 7 on the boundary, our results are in 
qualitative agreement with those presented in Ref. [pT]| . 
Comparison has also been made between perturbative 
results p^ , the numerical computations described here 
and the experimental data Q . The magnitude of 7 used 
in the perturbative analysis is the prime factor causing 
the mismatch between the perturbative and the numer- 
ical results. When 7 is large (=0.5 for our case), the 
first order perturbative approach is not expected to be 
accurate. 

While the results of the perturbative analysis and the 
numerical study are different quantitatively, they possess 
similar qualitative features, namely the onset of the ex- 
cluded angle ^'o as the domain size Rq increases, and 
then reaches a maximum of and then decrease as Rq 
continues to increase. These match the qualitative fea- 
tures that are present the experimental data shown in 
Fig. ||(b) . Experimental results is not reproduced in the 
numerical calculations when Rq is small. The discrep- 
ancies between the experimental data and the numeri- 
cal result imply that other interactions, neglected in the 
model, may be significant. 

We have also extended the numerical algorithm to the 
problem of bubbles. It is found that the transformation 
r — > r' = 1/r results in a new domain problem which al- 
lows us to solve the equilibrium conditions for the bubbles 
at the same level of efficiency and accuracy as those for 
the domains. Not only have we obtained results that are 
consistent with those in the perturbative analysis iQ,^ , 
we have also analyzed the effect of elastic anisotropy on 
the textures and the boundary of the bubbles, a task 
that is algebraically formidable in the perturbative anal- 
ysis. The influence of elastic anisotropy on the bound- 
ary is small while it significantly modifies the textures. 
This provides a means for the qualitative determination 
of the elastic anisotropy by observing the texture around 
the bubbles. The agreement between the perturbative 
results and the numerical computation is excellent. This 
is not surprising as a2 is not involved. The Rq depen- 
dence of is similar to the case of the domain, except 
that the maximum of ^0 is much smaller. The perturba- 
tive result agrees reasonably with the experimental data 
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as reported in Ref. 0. The numerical results match as 
well. 

In conclusion, we have successfully implemented a nu- 
merical algorithm that e nabl e us to analyze unambigu- 
ously a simple model, Eq.(2.1), of tilted ordered media in 
a non-trivial geometry imposed by experimental observa- 
tions. Using this numerical algorithm and its extensions, 
we are able to address the long-standing debate with re- 
gard to the origin on the cusp-like features observed in 
domains of Langmuir monolayers using an elastic model. 
Within the context of this simple model which addresses 
only the competition between the bulk elastic energy and 
the boundary energy, many qualitative features of the ex- 
perimental observation have been captured. Discrepan- 
cies cannot be avoided, as the real system is much more 
complex. The model we adopted has neglected other ef- 
fects and interactions that are present in real system, 
such as dipolar interactions and adjustments in the tilt 
degree of freedom. A combination of these effects may 
account for the discrepancies between the experimental 
data and the theoretical results. The apparently general 
numerical algorithm is, however, not capable of handling 
situation in which + ct" = at some points on the 
boundary. A different approach, such as the Wulff con- 
struction is required. Nevertheless, our numerical 
algorithm is versatile and can be extended to systems 
containing topological defects, or with the ordered phase 
filled in a nonsimply connected space. 
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APPENDIX A: VARIATIONAL FORMULATION 
OF THE FEM 

In fini te element analysis, we approximate Q(x, y ) by 
Eq. (2/7). The energy functional [6] given in Eq. (2.1) 
now becomes 



Hi@)^^ Jj\yQ\' + b[{-el + el) cos 28- 

29^85, sin 26]} dA <j> cr(i9 - 6) ds (Al) 

where we denote = (61, 82, ■ ■ ■ , 8jv)-^, and 8 = 
^(/?i8i. We differentiate Eq. (Al) with respect to a 
Qi yields 

— — — ^ K [ip,x(pjxil - 6 COS 28) + 

LPiyiPjy{l -I- &COS28) -1- 



h{(pix^jy + Viy^jx) sin26]c!^ej 
-Kb [ [(-8^ + 8^) sin 28 + 

28j;8j^cos28].^jdA - j) a'{d - e)ip,ds (A2) 

The equilibrium condition gi ves A Q = b with A and b 
provided in Eqs. (2.1C) and ( ^.11 ). 



APPENDIX B: INTEGRATION OVER A 
TRIANGULATION 



The integrals in Eqs. (2.10) and (2.11) over fl are 
broken up into sums of integration over the triangles in 
the triangulation of f2. Integration over the interior in- 
dividual triangle can usually be carried out analytically 
depending on the specific forms of the basis functions 
ipi{x,y) and the matrix elements Aij and bi. We have 
chosen ipi{x, y) to be a continuous, piecewise linear func- 
tion in X and y within a triangle. The line integral ds 
in bi will must be evaluated numerically, because inte- 
grand depends on the polar angle Lp, which is not lin- 
ear in X or y. This does not degrade the efficiency of 
the computation because first of all, only triangles whose 
perimeters coincide with T contribute to the line integral 
and secon dly it is a line integral over a short distance. 

In Eq. (^.7|), we express a function f{x,y) for {x,y) S 
n in terms of its values at the nodes of the triangulation 
and the corresponding basis functions ipi{x,y). Within 
an individual triangle we can write 



(BI) 



i=l 



where Ki is the index of the ith vertex of the triangle 
K, fxi = f{xK,,yKi), {xKi.yKi) are respectively the 
fimctional value of /(a;, y) and the coordinates of the ith 

vertex. ip'^\x,y) is the restriction of ipKi(x,y) in K. 
The actual index of the ith vertex is Ki. It is, however, 
awkward to carry the K in the symbol Ki throughout 
the discussion. We will use i to identity the vertex for 
simplicity from now on, i.e. is simplified as fi. We 
introduce a set of natural coordinates u and v such that 



(B2) 



where u e [0, 1], w € [0, 1] and u + v <\. Transformation 
between variable sets x — y and u — v can be obtained 



from Eq. (B2) by substituting / with x and y. We then 
have the foUowings relations 

" = ^ [(2/3 -yi)x- (x3 -xi)y + xsyi - xiys] , (B3) 
'K ^'■^^ ~ 2^2) a; - {xi -X2)y + xiy2 - X2yi] , (B4) 
where A is the Jacobian determinant given by 
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A = 



d{x,y) 
d(u, v) 



Identifying Eqs. (Bl) and (B2), we find ipi{u,v) — 1 
have 



xiy2 - yiX2 + X2y3 - X3y2 + x^yi - xiy^. (B.5) 
— f , ip2iu, v) — u and 933(1*, v) — v. In terms of x — we 



fiix,y) = 


1 

A 


[{y2 




{X2 


- X3) y + 0:2^3 


- 2:32/2] , 


'P2{x,y) = 


1 

A 




~yi)x- 


{X3 


- xi)y + x^yi 


- xm] , 


f3{x,y) = 


1 

A 


[(yi 


-2/2)2;- 


{xi 


- X2)y + Xiy2 


- 2:22/1] ■ 



(B.6) 
(B.7) 
(B.8) 

Evaluation of the matrix element Aij involves of the following area integrals which can be computed analytically. 
The trivial one is the area of K, which is J^, dA = |A|/2 and 



cos 2Q dA 



K 



sin 29 dA 



K 



4 
4 



cos28i 



cos 202 



cos 283 



(61 



e2)(e2 

sin 261 



63 



(62 



e3)(e3 

sin 262 



61) (63 



ipi cos 20 dA = 



K 



•^2 COS 28 dA — — — 



K 



. (61 - e2)(e2 - 63) (62 - e3)(e3 - Oi) 

as Lpix and ipiy are constants. Evaluation of hi involves 

T 
T 

K ~ 



ei)(ei-e2) 

sin 263 

(e3-ei)(ei-e2) 



ip3 COS 20 dA 



K 



ipi sin 20 dA 



K 



ip2 sin 20 dA 



K 



(^3 sin 20 dA 





COS 201 




sin 201 — sin 203 


sin 201 — sin 202 


(63 


-0i)(0i- 

cos 202 


02) 


^ 2(01-03)2(03-02) 

sin 202 — sin 203 


2(01-02)2(03-02) 

sin 202 — sin 201 


(01 


-02)(02- 

cos203 


63) 


+ 2(02-03)2(03-01) 

sin 203 — sin 202 


2(01 -02)2(03-01) 

sin 203 - sin 201 


(02 


-e3)(03- 

sin20i 


0i) 


2(02-03)2(02-01) 

COS 201 — COS 203 


2(01 -03)2(03 -61) 
cos 201 — cos 202 


(03 


-0i)(0i- 
sin 202 


02) 


2(01-03)2(03-02) 
COS 202 — COS 203 


2(01 -02)2(03-02) 

cos 202 — cos 201 


(01 


-02)(02- 
sin 203 


63) 


2(02-03)2(03-01) 
COS 203 — COS 202 


+ 2(01-02)2(03-01) 

cos 203 — COS 201 


(02 


-03)(03- 


0i) 


2(02 -03)2(02-01) 


+ 2(01 -03)2(03-01) 



(B.9) 
(B.IO) 

(B.ll) 
(B.12) 
(B.13) 
(B.14) 
(B.15) 
(B.16) 



The above formulae will work only if 0i 7^ ©2 7^ 03- Let us consider cases where the values of = 0^; at 2 vertices 
of triangle K and = 0o at the other vertex. We obtain the following for the integrals in A^ 



K 



K 



cos 20 rfA = — - 
4 

4 



sin 20 dA 



COS20£; — cos20o 

(0iJ - 0O)2 

sin 20^; — sin 20o 

(0B-0O)2 



2sin20£; 



0B-0O 

2COS20£;' 
0B-0O. 



(B.17) 
(B.18) 



We denote the restrictions of the basis functions at the nodes that have = 0^;, and ^po the restriction of the 
basis function at the node that has — ©q. We arrive at 



K 



K 



|A| 

(i3£;COs20 dA= — 

4 

|A| 

ipo cos 20 = 



sin20£; 



COS20£; 



l20b-: 



L20f 



LPE sin 20 dA ~ 

4 



Lpo sin 20 dA 



K 



4 



0B - 0O (©£ - ©o)2 2(0B - 0o 

cos20B + cos20o sin20£; — sin20o 

(0s-0o)2 ^ {Qe-^o? 
COS 2Qe sin 20 e cos 20^ — cos 20o 

"0B-0O ^ (0B-0o)2 ^ 2(0B-0o)3 

sin 20£; + sin 20O cos 20^; — cos 20o 



(0i;-0o)2 (0B-0o)3 

for the integrals required to evaluate hi. Finally, when 0^ = 05 for all i's, one will need ^pE dA — |A|/6. 



(B.19) 
(B.20) 
(B.21) 
(B.22) 



12 



APPENDIX C: DERIVATIVES ON THE 
BOUNDARY 

One of the biggest benefits of the FEM is that it en- 
ables straightforward determination of the derivatives of 
6 on the boundary. The tangential derivative at node i 



dt 



1 



2e''Vl + 



(C.l) 



The normal derivative of Q is given by Eq. (2.5) which 
reads 



dQ _ e^a'i'd - 9) + nbQt sin 2(6 - d) 
'dn ~ K[l-6cos2(e-i?)] ' 



(C.2) 



[1] CM. Knobler and R.C. Desai, Annu. Rev. Phys. Chem. 
43, 207 (1992). 

[2] J. Ruiz-Garcia, X. Qiu, M.-W. Tsao, G. Marshall, CM. 
Knobler, G. A. Overbeck and D. Mobius, J. Phys. Chem 
97, 6955 (1993); D.K. Schwartz, J. Ruiz-Garcia, X. Qiu, 
J.V. Seligner and CM. Knobler, Physica A 204, 606 
(1994);J.V. Selinger, Z.-G. Wang, R.F. Bruinsma and 
CM. Knobler, Phys. Rev. Lett. 70, 1139 (1993). 

[3] X. Qiu, J. Ruiz-Garcia, K. J. Stine, CM. Knobler, J.V. 
Selinger, Phys. Rev. Lett. 67 703 (1991). 



[4; 
[6; 

[7 
[8 
[9 
[10 

[11 

[12 
[13 

[14 

[15 
[16 
[17 

[18 



T.M. Fischer, R.F. Bruinsma and CM. Knobler, Phys. 
Rev. E 50, 413 (1994). 

S. Riviere and J. Meunier, Phys. Rev. Lett. 74, 2495 
(1995). 

N.D. Mermin, in Quantum Fluids and Solids, edited by 
S.B. Trickey, E. Adams and J. Duffy (Plenum, New York, 
1977). 

J. Fang, E. Teer, CM. Knobler, K.-K. Loh and J. Rud- 
nick, Phys. Rev. E 56, 1859 (1997). 

S.A. Langer and J.P. Sethna, Phys. Rev. A 34, 5035 
(1986). 

I. Kraus and R.B. Meyer, Phys. Rev. Letts 81, 3815 
(1999). 

J. Rudnick and R. Bruinsma, Phys. Rev. Lett. 74, 2491 
(1995). 

P. Galatola and J.B. Fournier, Phys. Rev. Lett. 75, 3297 
(1995). 

J. Rudnick and K.-K. Loh, Phys. Rev. E 60, 3045 (1999). 
C Johnson, Numerical solution of partial differential 
equations by the finite element method (Cambridge Uni- 
versity Press, Cambridge, 1987). 

W.H. Press, B.P. Flannery, S.A. Tuekolsky and W.T. 
Vetterling, Numerical Recipes in C, The Art of Sci- 
entific Computing (Cambridge University Press, Cam- 
bridge, 1991). 

K.-K. Loh and J. Rudnick, Phys. Rev. Lett. 81, 4935 

(1998) . 

D. Pettey and T.C Lubensky Phys. Rev. E 59, 1834 

(1999) . 

E. Teer and CM. Knobler, private communications; C. 
Lautz, T.M. Fischer, M. Weygand, M. Losche, P.B. 
Howes and K. Kjaer, J. Chem. Phys. 108, 4640 (1998). 
J. Fang and CM. Knobler, unpublished. 



13 



